# ---------------------------------------------------------
# TUM - Technichal University of Munich
#
# Original Authors: Taylor Jones
#                   Magdalena Altmann
# Further Authors:  Adrian
# Date: 2020-05-06
# Purpose: This code is only to plot the BIMs. It has been
#          simply copied from the existing Inversion code.
# ---------------------------------------------------------

rm(list = ls(all.names = TRUE))  # clear R environment (clear workspace)
graphics.off()  # clear plots
# dev.off()  # clear plots

library(reshape2)
require('raster')
require('leaflet')
require('proj4')
require('ncdf4')
require('xts')
require('rasterVis')
require('gridExtra')
require('ggplot2')
require('viridisLite')
require('MCMCpack')
require('Metrics')
require('rsq')
require('base64')
require('base64enc')

t1 <- Sys.time()
#dates <- c(20180818,20180819,20180820,20180821,20180822,20180827)

dates <- c(20180827)
#dates <- c(20180818, 20180820)
#dates <- c(20180819, 20180821, 20180827)

ems     <- c('mb','ma','ki','kj','wa') #'mb','ma','ki','kj','wa'

#ems     <- c('mb','ki','kj','wa') 
#ems     <- c('ki')

# get directory of source file and set it as working directory
path_inverse_modeling_engine <- dirname(rstudioapi::getSourceEditorContext()$path)  # determine source file path
setwd(path_inverse_modeling_engine)  # set path as working directory


source('plotting.R')             # has all the plotting functions
source('load_inventories.R')     # has functions to load different inventories into raster stacks
source('inversion_functions.R')  # The rest of the functions

#Constants and Conversions:
#to_ggyr <- 1e-6 #16.04*365.25*60*60*24*1e-9  #umol/m2s -> Gg/yr ...still need to multiply by area in km2!
to_ggyr <- 16.04*365.25*60*60*24*1e-9  #umol/m2s -> Gg/yr ...still need to multiply by total domain area in km2!

# using preprocessed data: apply_scaling = FALSE; using observations directly: TRUE
apply_scaling = TRUE
# solve inversion for all measurement days together: invert_all_days = TRUE
invert_all_days = FALSE


#Parameters for Inversion:  -------------------------------------------------------------------------------------------------------------
sigma_sector_priors <- list("total"=1e2)
sigma_obs_prior     <- 2e-4      # ppm  (0.2 ppb)
sigma_bkgd_prior    <- 1e-3      # ppm  (10 ppb) 1e-2 
t_back              <- 3*60*60   # seconds: The timescale (e-folding rate) of background variability
bkgd_prior          <- 1.85      # ppm: prior background value

#Parameters for Bootstrap: --------------------------------------------------------------------------------------------------------------
sector_to_bootstrap <- "total"
bs_x_err            <- 2e-4  # ppm:  (0.2 ppb) error with which to draw for bootstrap
bs_y_err            <- 2e-4  # ppm:  (0.2 ppb) error with which to draw for bootstrap
bootstrap_reps      <- 5000  # times to run each of the bootstraps

#Load the Prior: ------------------------------------------------------------------------------------------------------------------------
date         <- dates[1]
foot_file    <- paste0("foot/",prefix,"_",date,"_ERA5_bw8_small.nc") # footprint files (needed here just for the projection)
loaded_prior <- load_inv(foot_file)                                  # read the prior maps and reproject it to the footprint projection.
sectors      <- loaded_prior$sectors                                  # list out the IPCC sectors listed in those maps
ras_stack    <- loaded_prior$ras  
prior_totals <- cellStats(ras_stack*(raster::area(ras_stack)),"sum")*to_ggyr    # Add up all the emissions and convert to GG/yr
names(prior_totals) <- sectors                                        # Name these totals with the sector names


#I'll need these later:
out_df_total <- data.frame()
bk_df_total  <- data.frame()
bars_df      <- data.frame()
y_df_total   <- data.frame()
bim_df_total <- data.frame()

bars_total   <- data.frame()

r2_enhancement_total   <- data.frame()
r2_concentration_total <- data.frame()

x_hat_total            <- data.frame()


bim_df <- load_bims(foot_file,ems)


plt_BIM(bim_df, pltlog=T)

t2 <- Sys.time()
t_elapsed <- t2 - t1
message("Time elapsed: ", round(as.numeric(t_elapsed, unit="secs"), digits=2), " s")
